A tangential force model for interactions between bonded colloidal particles 
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Recently, Pantina and Furst (Phys. Rev. Lett., 94(13), 138301, 2005) experimentally demon- 
strated the existence of tangential forces between bonded colloidal particles and the capability of 
these bonds to supporting bending moments. We introduce a model to be used in computer sim- 
ulations that describes these tangential interactions. We show how the model parameters can be 
determined from experimental data. Simulations using the model are in agreement to the measure- 
ment by Pantina and Furst. Application of the model to an aggregate with fractal structure leads 
to more realistic behavior than using classical approaches only. 

PACS numbers: 61.43.Hv, 82.70-y, 47.57.-s, 52.65.Yy, 87.15. A-, 



I. INTRODUCTION 



The structure of colloidal aggregates is important in 
various applications (e.g. pharmaceuticals, food process- 
ing, nano-particle synthesis). To address structural as- 
pects micro-simulation of aggregating colloidal particles 
are an important and growing field in colloid science. 
Micro-simulation of aggregates allows to investigate the 
structural evolution of aggregates by tracking the trajec- 
tories of the constituent primary particles. These tra- 
jectories are obtained from solving Newton's equation 
of motion for all the aggregates primary particles. In 
the future, insight in structure formation may be ex- 
ploited to tailor aggregate structures by optimal process 
design. In the literature some work can be found on sim- 
ulating aggregate structure evolution in hydrodynamic 
environments. Generally, one must distinguish the hy- 
drodynamic and the particle interaction forces. Initial 
attempts of structural modeling have been carried out 
by Doi and Chen [l|, Hj. For the hydrodynamic forces 
they used the so-called free draining approximation ac- 
cording to which the hydrodynamics forces on the par- 
ticles are strongly simplified. Each particle experiences 
only Stoke's drag force. Thus, all flow field perturba- 
tions induced by the particles, which potentially affect 
the flow around neighboring particles, are neglected. For 
the particle interaction they also used a simple model 
where sticking particles can roll on each other and the 
bond between two particles breaks if the normal forces 
exceed a critical value. Higashitani et al. [|[ performed 
simulations, where the hydrodynamic and inter-particle 
forces are considered in much more detail. Inter-particle 
forces were obtained by the classical Dcrjaguin, Landau, 
Verwey, Oberbeek (DLVO) theory For particles in 
contact they used the model of Cundall and Strack [5[ 
which is widely used in Discrete Element Modeling of 
granular matter [fjj] . 
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Regarding the hydrodynamic forces they accounted for 
screening of inner particles from the flow field by means of 
detailed geometrical computations. Similar studies have 
been done by Fanelli et al. 0, Hi who also used a Dis- 
crete Element Method (DEM) and DLVO forces to simu- 
late dispersions of colloidal aggregates. Harada et. al. @ 
examined the structural change of non- fractal clusters. 
To compute the hydrodynamic forces they used Stoke- 
sian dynamics [IcL [III ] which allow computation of the 
full, hydrodynamic interaction on the basis of the Stokes 
equation. As direct particle interactions they considered 
a retarded attractive van der Waals potential. Most of 
these studies assume normal forces between particles. 

There are only a few simulation studies where non- 
normal forces are included. In the work by Higashitani 
et al. || the contact model of Cundall and Strack [j| is 
employed. That model is designed to capture sticking 
and sliding friction. However, as it will be shown, the 
classical Cundall-Strack model is not capable of quali- 
tatively describing the experimentally observed behavior 
of bonded colloidal particles. In the field of disordered 
networks, e. g. gels and glassy structures, some models 
for tangential forces capable to supporting bending mo- 
ments were derived (see e. g. [t3. Il3l H^.Tl5j). Potanin 
[l|| adopted the main ideas of these models to apply 
them in the context of colloidal aggregates. He used a 
Hamiltonian used in Born's Model for the elasticity of 
microscopic networks fl7^ . This model was applied to 
the simulation of aggregates under static conditions but 
is not suited to describe the dynamic behavior, e. g. the 
rotation of an aggregate, correctly. In subsequent work 
[l|| , the author calculated the tangential forces based on 
a Hamiltonian derived by Kantor and Webman [Ti l] . This 
Hamiltonian based on a three body approach. The bend- 
ing energy is proportional to the variation of the angle 
between two neighboring bonds. With this approach the 
breakage behavior of aggregates formed by diffusion lim- 
ited cluster aggregation (DLCA) was investigated. How- 
ever, this approach does not allow irreversible rearrange- 
ment of particles. Thus, the approach is restricted to sit- 
uations where restructuring effects are irrelevant. West, 
Melrose and Bell [lj| presented another approach to in- 
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elude tangential and bending forces. They replaced the 
particles by a stiff trimer of particles and the basic col- 
loidal bond was taken as a superimposition of 3x3 sphere 
interactions. Botet and Cabane [2(| introduced a model 
where the bond between two spherical colloidal particles 
is described by springs which connect pins on the spheres 
surface. These pins are randomly distributed over the 
spheres. Whenever the distance between two pins on dif- 
ferent spheres is smaller than a characteristic threshold, a 
spring between these pins will be initiated. The two pins 
cannot form another bond as long as the bond is present. 
If a spring exceeds a maximum elongation the spring will 
be destroyed. That springs causes both normal and tan- 
gential interaction between the spheres. A comparison of 
that model to the one proposed in this work will be given 
in section ITV CI 

Recently, experimental evidence of bending moments 
has been published by Pantina and Furst [2l| . They mea- 
sured that even a single bond is capable of supporting a 
bending moment. Except the model of Botet and Cabane 
[20l |. none of the above mentioned models can describe 
both observed effects. In this paper we present a model 
to be used in DEM simulations which is able to describe 
the phenomena found by Pantina and Furst [2l| . We 
furthermore show that all necessary model parameters 
can directly be obtained from their experiments. We will 
show that in some special cases our model reproduce the 
Hamiltonian used in [l4| and we will compare our model 
to the models mentioned above. 

The organization of the paper is as follows. Section 
HI! summarizes the basic experimental findings of Pantina 
and Furst regarding tangential forces as they are of major 
importance to the model developed in this contribution. 
In section IIIII we will briefly revisit the classical DLVO 
forces. Section HV1 introduces the novel model for tangen- 
tial forces, presents the method to determine the model 
parameters and compares our model to previous ones. 
Section [V] explains the basic simulation technique and 
shows simulation results. Finally, section IVT1 summarizes 
our findings and draws some conclusions. 

II. EXPERIMENTAL OBSERVATIONS BY 
PANTINA AND FURST 

In their experiments Pantina and Furst [2l| used a lin- 
ear chain of polymethylmethacrylate (PMMA) particles 
immersed in an MgCl2 aqueous solution. The terminal 
particles of the chain were fixed by optical tweezers and 
the mid particle was pulled by another optical tweez- 
ers perpendicular to the chain direction. If only central 
forces acted between the particles the formation of a tri- 
angular structure can be expected. However, it turned 
out that the positions of the particles are in good agree- 
ment with the shape expected from a thin elastic rod: 

y( x ) _ 1 ( L x 2 M 3 \ m 
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Here y(x) is the deflection as function of the position x, 
L is the length of the aggregate, E is the Young modulus 
and I the second moment of area. This is clear evidence 
that bonds between colloidal particles are capable of sup- 
porting bending moments. Furthermore, they measured 
the bending rigidity, k, defined as the constant of propor- 
tionality between the deflection 6 of the aggregate and 
the applied force: 

Fbcnd = kS . (2) 

Additionally, they found that n oc L~ 3 as expected from 
equation (TT]). The bending rigidity can be expressed by 
the relation 

(d \ 3 
LJ ' ^ 

where a is the particle radius and kq is the bending 
rigidity per bond. They measured these quantities as 
functions of the MgCl2 concentration. Furthermore they 
found that there is a critical bending moment M c , above 
which the particle begin to slide and rearrangements of 
particles occurs. 

Pantina and Furst [2l[ presented a possible explanation 
for the tangential forces in terms of the Johnson Kendall 
Roberts (JKR) theory for adhesive surfaces. They re- 
lated the single bond rigidity kq to the work of adhesion 
W s i derived from the JKR theory and to the Young mod- 
ulus of the particles. In subsequent work, Pantina and 
Furst j22} generalized that theory to the case that diva- 
lent ions are present. In that case ion bridges between the 
spheres appear which increases the attraction between 
the particles and leads to a higher bending rigidity. 

III. DLVO FORCES 

Let us briefly revisit the DLVO theory. The first ingre- 
dient of the DLVO theory are Van der Waals forces. In 
the framework of the non-retarded Hamaker approxima- 
tion, the mutual interaction potential between two par- 
ticles can be found in standard textbooks (e. g. pi 1231]): 

T „_, A f 2a 2 2a 2 , f R 2 - Aa 2 \\ 

where R is the center-center distance between two parti- 
cles and A is the Hamaker constant, depending on parti- 
cles' and fluid's properties. The second ingredient is the 
electrostatic double layer theory. Here we use the Der- 
jaguin approximation with the assumption of constant 
surface potential. The mutual interaction potential can 
again be found in the literature [J, [23| : 

V(R) = ln[exp(-A(i? + a))] . (5) 

Here e is the electric permittivity of the carrier fluid, 
</>o is the surface Potential and A is the Debye-Hueckel 
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parameter defined as 
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where e is the elementary charge, /cb is the Boltzmann 
constant, T is the temperature, rij is the ion concen- 
tration of the i'th ion species and Zi is the correspond- 
ing valency. The inverse of the Debye-Hueckel param- 
eter is a measure for the magnitude of the screening 
length for electric fields in an electrolyte solution. Be- 
sides these standard ingredients we use an additionally 
repulsive short range force (Born repulsion force), mak- 
ing sure that particles cannot collide. We use a formula 
derived by Feke et al. [24| : 
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FIG. 1: Interaction potential between two equal spherical 
particles of radius a — 0.735 ^m. The potential comprises 
attractive Van der Waals interactions with a Hamaker con- 
stant of A = 0.062 eV, repulsive electrostatic interactions in 
an 150 mM MgCl2 aqueous solution, a surface potential of 
(7) c/>o = 40mV and the Born repulsion where N is assumed to 
be N = 10" 23 



In this expression R = R/a is the ratio of the center- 
center distance and the particle radius. As discussed by 
Feke et al. [H| N lies in the interval lCT 18 to 1CT 23 . In 
our simulations we used N — 1CP 23 . Figure Q] shows the 
interaction potential between two particles for conditions 
corresponding to the experiments by Pantina and Furst 
pl| . As expected for this parameter set, the aggregation 
of particles is not hindered by an energy barrier and the 
colloid is in the regime of fast coagulation 0] ■ The equi- 
librium surface-surface distance, that is the position of 
the potential minimum, is typically in the order of some 
Angstrom. 



IV. NOVEL TANGENTIAL FORCE MODEL 
A. The model 

There is quite a number of tangential force models 
available for DEM simulations. For example the mod- 
els by Haff and Werner (25|, Walton and Brown [26| or 
Cundall and Strack [f| . For some detail on these models 
the reader is referred to 27]. For this work the main de- 
ficiency of all these models is that most of them are not 
capable of supporting bending moments between bonded 
particles, except the ones which will be discussed and 
compared to our model in section HV CI In the following 
we propose a novel tangential force model, similar to that 
by Cundall and Strack, which can reproduce the phenom- 
ena described in section [TTJ In the Cundall-Strack model 
a spring £y with rigidity k t will be initialized when two 
particles, i and j, get into contact. This spring grows pro- 
portional to the relative tangential velocity at the contact 



point: 
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t = v t 



(8) 
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Here to is the time when the particles get into contact. 
The relative tangential velocity is given by 



v t = (vj - Vi) t + a (u>j + u>i) x rii 



(9) 



where n,j is defined as (rj — ri)/\rj — r,-|. The subscript 
t denotes the projection of the relative velocity onto the 
tangential plane. If the force due to the tangential spring 
exceeds an upper bound, which in the work by Cundall 
and Strack Q is the Coulomb friction fiF n , the spring 
will be set to 



(10) 



where [i is the friction coefficient. To make sure that 
the tangential force acts only in tangential direction, the 
Cundall-Strack spring is mapped to the tangential plane 
perpendicular to after each time-step. The tangential 
forces and torques acting on the i'th and the j'th particle 
are given by 



Mi = RiF t A x n t 



-RjFt,j x m. 
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This model is able to capture a lot of phenomena like slid- 
ing friction and sticking friction. However, it is not able 
to support a bending moment between two bonded par- 
ticles, which becomes obvious by the following example. 



FIG. 2: Example for the bending torque, the lower particle 
is fixed and on the upper particle acts a tangential force (left 
side). If the bond between the particles is able to support 
a bending moment, a stationary angel (j> < n/2 between the 
original contact line and the stationary contact line should be 
reached. 



FIG. 3: Sketch of tangential force model, (a) When two parti- 
cles get in contact two thought rods are initialized. The rods 
are rigidly connected to the center of one particle and reach 
the other particles center, (b) According to the tangential 
movement of the particles the springs and will elon- 
gated. The arrows denote the direction of the vectors £^ and 



We assume a fixed (no translation or rotation) particle 
bonded to a second particle as sketched in figure An 
external force perpendicular to the center-center line of 
the particles acts on the second particle. Now we look 
for a static solution of the evolution equations, in which 
all resulting forces and moments have vanished. From 
equations (fTTj) and (fT2|) immediately follows that this is 
only the case if £ — 0. This in turn means that the whole 
external force must be equilibrated by the normal forces 
between the particles and this can only be the case if 
the center-center line is finally parallel to the direction 
of the external force. However, this means that the bond 
cannot support any bending moment. In order to derive 
a similarly simple model, however capable of supporting 
bending moments, we use the following reasoning: When 
two particles get in contact, we introduce two thought 
rods rigidly connected to one particles center and reach- 
ing to the center of the other particle as sketched in figure 
[31 Between the end point of the rod and the center of 
the other particle a spring will be initialized. This spring 
grows proportional to the relative tangential velocity be- 
tween the rods end points and the particle center. The 
evolution equations for the springs are then 



k%j = ( v j ~ v i)t ~ (w< x mj) 2a , (13) 
kji = (vi - Vj) t + (uj x mj) 2a , (14) 

where £^ ■ , £ ^ are defined in figure [3] The forces and 
moments acting on the particles are therefore 

Fi = h - Cji) , Mi = 2ak t n ij x £ tf , (15) 
Fj = h (£., ~ £.,) , M 3 = -2ak tnij x t,„ • (16) 

Equivalent to the Cundall-Strack model both tangential 
springs are mapped to the plane perpendicular to fly 
after each time step. The springs stop extending if its 



elongation exceeds a maximum value £ max - Unlike the 
model of Cundall and Strack, this model is able to sup- 
port bending moments. Let us demonstrate that by the 
same example as above. 

By solving the steady state equations one find in the 
case of small deflections for the reorientation angle </>: 



k f 2a 



(17) 



Note that our model has some similarity to a theoreti- 
cal approach introduced in [28[ where a random network 
of springs (normal and tangential) was used to predict 
the elastic moduli for disordered packings of intercon- 
nected spheres. Here the tangential spring stiffness was 
determined by the predictive model of Pantina and Furst 

mm 

We assume the model to be valid also for non-tenuous 
aggregates as long as pairwise contact forces can be 
assumed. In this contribuition we content ourselves 
with polymeric primary particles as investigated in [2lj . 
Whether our model is also applicable to other types of 
particles (e.g. inorganic and surfactant-coated colloids) 
can not be revealed by the model itself but has to be 
clarified experimentally. 



B. Parameter determination 

The introduced model contains two parameters. The 
spring stiffness k t and the maximum spring length £ max . 
Both parameters can be determined by the experiments 
presented in [2l[ . As shown in appendix [B] the static 
shape of a linear chain of particles under a bending stress 
is given in leading order approximation by 
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where L is the center-center distance of the first and the 
last particle in the chain. By comparison with equation 
(TTJ one finds 



EI 



(19) 



From the elongation of the chain, 8 = y(0) — y(L/2), one 
finds by comparison with equation @ K = 192(a/ L) 3 kt 
and therefore from ([3]): 



h = 



HQ 

192 



(20) 



Using equation (|20p one directly obtains the model pa- 
rameter for the tangential stiffness k t from the measured 
value kq. The value of £ max can be obtained from the 
measurement of the critical bending moment M c in [2l| . 
The bending torque acting on a particle is given by equa- 
tion ([T5]l or ([T5]), Since n,j and £ i3 are perpendicular to 
each other £ max has to be 



Mc_ 
2akt 



(21) 



Thus, all parameters of the force model can be deter- 
mined from experimental data. Alternatively, the pre- 
dictive formulas derived by Pantina and Furst [2ll |22| 
for the single bend rigidity kq may be used to estimate 
the model parameters. 



C. Comparison to other models 

In [l8| a model derived by Kantor and Webman [l4| 
was applied to computer simulations of two dimensional 
colloidal aggregates in shear flows. That model assumes 
that the contribution of tangential deformation to the 
strain energy is proportianal to 



E oc J2 . 



(22) 



where 8<&ijk is the change in the angle between the bonds 
and (i,fc) connected to particle i. The contribution 
to the interaction force can be described by 



ijk^ik - ( n ij - n i. 



(23) 



It is mentioned in [l8| that a chain of particles connected 
by such forces behaves like a thin elastic rod with the 
same length. As shown in appendix [3] our model leads 
in the case of a linear structure (i. e. every particle is 
bonded to a maximum of two other particles) to an elastic 
equilibrium energy of 



E = a 2 k t s ®i 



(24) 



where <5$i denotes the change of the angle between the 
two bonds ending on particle i. Due to the fact that 



and (|22j) are equivalent, the equilibrium elastic deforma- 
tion of our model is also equivalent. 

The main difference between the two models is that 
our model takes only pairwise interactions into account, 
while the model by Kantor and Webman [l4[ is a three 
particle interaction model. The tangential forces and the 
bending energy are there related to the change of angles 
between neighboring bonds. However, one would expect 
that the energy is stored directly in the bonds. It remains 
unclear why the angle between two different bonds is the 
measure for the elastic energy. Although the same be- 
havior is observed in our proposed model, the reason is 
quite obvious. The energy is directly stored in the bonds 
but neighboring bonds are related due to the equilib- 
rium conditions. In contrast to the model of Kantor and 
Webman the proposed model is suited to describe non- 
elastic displacements of the bonds, which occur if the 
maximum supported bending moment is reached. Fur- 
thermore, the tracking of angles between all bonds ending 
on the same particle may cause much higher implemen- 
tatory (and possibly computational) effort compared to 
the proposed model. 

To avoiding the computational problems caused by 
tracking of the angles, West, Melrose and Ball [l!| pro- 
posed a model where the basic colloidal unit is replaced 
by a trimer of stiff spheres and the interactions are as- 
sumed as superposition of 3 x 3 pairwise interacting cen- 
tral forces. Depending on the relative orientation of the 
trimers, bonds between two trimers are able to support 
bonding moments. Contrary to our proposed model this 
only provides qualitative information. Besides, no rela- 
tion between experimental data and model parameters 
is available. Furthermore replacing one particle by three 
particles leads to a higher computational effort. Never- 
theless, the ansatz of West, Melrose and Ball is probably 
superior to our model for the simulation of colloids con- 
taining non-sphcrical primary particles. 

The model of of Botet and Cabane [2(| is in princi- 
ple able to capture the effect observed by Pantina and 
Furst 21]. They modeled the interactions between two 
spheres by a number of springs connected to pins which 
are randomly distributed over the spheres. In compar- 
ison to our model the model of Botet and Cabane [2(| 
consumes much more computational resources in order 
to handle the pins and springs. Furthermore, the estima- 
tion of the model parameters, i. e. number of pins, spring 
stiffness, equilibrium and breakage lengths, seemed to be 
more difficult. 



V. SIMULATION METHOD AND RESULTS 

We assume spherical mono-disperse colloidal particles 
immersed in an aqueous MgCl2 solution. Unless other- 
wise noted we used the parameters from [2l|. The ra- 
dius of the particles is a = 0.735 fim. The surface po- 
tential is 4>o = 40 mV. We use the Discrete Element 
Method 0, @ to simulate aggregate structure evolution. 
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The main idea of this method is to track the trajectories 
of all primary particles by solving Newton's equations of 
motion numerically. In general the state of the particle 
system is given by the particle positions {ri, r 2 ...rjv}, 
velocities {v\, v%, v „} and by the angular velocities 
{u)\,u)2> "Wn}. Note that there is no need to track the 
particles' orientation angles as we content ourselves with 
spherical particles. For the interactions between the fluid 
and the particles we use the free draining approximation: 
Every particle experiences the unperturbated flow field as 
if no other particle were in the flow. The drag force and 
torque act ing on the «'th particle is then given by Stokes 
formulas: [29l ] 

Fdra g ,i = 6m]a(vi - ufo)) , (25) 
Mdrag.i = 87rrya 3 (o; t - n(n)) . (26) 

With r), Vi, and «(/%) being the dynamic viscosity of the 
carrier fluid, the velocity of i'th particle and the fluid 
velocity at the position of i'th particle, respectively. £1 = 
I/2V x u is the vorticity of the fluid velocity field. For 
the considered particles, effects of inertia can surely be 
neglected. Therefore, the particle dynamic is governed 
by the overdampded equations of motion: 
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Fi + uiri) . 
■Mi + n(n 



(27) 
(28) 



The forces Fi and the torques Mi include all interaction 
effects acting on the i'th particle. For solving the equa- 
tions of motion we use Heun's method, which has a global 
error of order At 2 and is similar to the velocity- Verlet 
method, often used in molecular dynamics and DEM sim- 
ulations [30| . Note that the velocity- Verlet method itself 
is explicitly designed for solving Newton's equations of 
motion and cannot be used in overdamped dynamics. In 
principle, we can use the force models described above to 
simulate the colloid. However, from the slope of the po- 
tential plot (figure [1} a computational problem becomes 
apparent. In the neighborhood of the potential minimum 
the interaction forces change rapidly on very small length 
scales. Therefore it is necessary to solve the equations of 
motion with very small time steps to track the details of 
motion and avoid instabilities of the numerical solution. 
We found that the magnitude of time-steps must approxi- 
mately be 10~ 9 s to track particle motion correctly. How- 
ever, if two particles are bonded to each other, the center- 
center distance remains approximately constant and only 
the angular orientation of the particles can change sig- 
nificantly. In order to avoid the need of such small time- 
steps we introduce a constraint if the distance between 
two particles i and j becomes smaller than a critical dis- 
tance d c . The equations of motion are then solved with 
the constraint |rj — rj\ = d c . For solving the constrained 
equations of motion, we adopted the RATTLE algorithm 
derived by Andersen (3fj , which fulfills the constraint up 
to a given tolerance tol, to Heun's method. In this paper 
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FIG. 4: Simulated shape of a bended Il-particle aggregate for 
different material parameters. The used tangential stiffness 
are kt = 0.785 mN/m (circles), kt = 1.1 mN/m (crosses), 
k t — 1.7 mN/m (squares) and kt = 3.4 mN/m (diamonds). 
The Symbols are simulation results and the lines are obtained 
from (Tj]). 



we used d c = 1.1 run and tol — 0.1 nm. By careful inves- 
tigation of numerical results we found that the value of 
d c has a negligible effect on the results to be presented 
below as long as d c is much smaller than the primary par- 
ticles' diameter a. With the chosen parameters we were 
able to use time-steps of the order of 10~ 6 s, which leads 
to remarkable improvements in simulation time. How- 
ever, the time needed by the Anderson algorithm is pro- 
portional to the number of bonds in the aggregate and 
especially for large compact structures the computation 
time may grow fast with the number of particles used in 
the simulation. 

In order to reproduce the results of Pantina and Furst, 
we performed simulations where a linear chain of par- 
ticles is bended. We applied an external force Feend 
directed perpendicular to the chain of particles to the 
middle particle. We compensate this force by applying 
— Feond/2 on the terminal particles of the chain. Then 
we run the simulation until the static shape is achieved. 
Note that the value of the fluid viscosity has no influence 
on the static solution, but only on the time needed to 
achieve it. Figure 2] shows simulated deflection curves for 
an aggregate comprising 11 particles for different spring 
stiffness k t — 0.69; 1.1; 1.7; 3.4 mN/m. These values 
correspond to the measured kq for different MgCl2 con- 
centrations (150; 250; 375; 500 mM, respectively 
All four curves are well described by equation ([1]) (solid 
lines) with EI obtained from equation (J2TJJ) . 

Figure [5] shows typical deflection vs. force curves. Here 
the measured value of n$ = 0.13 N/m for a 150 mM 
MgCl 2 concentration was used [2l|. For small deflections, 
the aggregate is bended similar to a rigid rod as described 
above. In that regime the deflection vs. force curve shows 
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FIG. 5: Typical deflection-force curves for a 11 particle (u 
line) and a 23 particle aggregate (lower curve). The pa 
eters for the tangential force model are k — 0.69 niN/m 
£max = 30.48 nm obtained from measured data in [211 ] : 
150 MgCb solution. The dashed line shows the experii 
tally obtained linear part of deflection-force curve taken 
a data plot (figure 2) presented by [2lj . 



a linear behavior. If the bending force exceeds a cri 
value, rearrangement of the particles occurs. Beyond 
point, the particles form a near triangular structure, 
simulation results in the linear regime and the cri 
force where the rearrangement occurs are in agreer 
with the experimental data presented in [2l[ figui 
Note that for experimental reasons the relevant pai 
the data presented by |2l| are on a deflection int< 
of approx. 1.5 — 2 fim. The reason is some tortou 
in the chains, which must be unbend before the linear 
elongation regime starts [35[. The simulated behavior 
after reaching rearrangement is different from the exper- 
imental data. This is due to the fact that the terminal 
particles in the experiment are trapped by optical tweez- 
ers while in the simulation only a force perpendicular to 
the linear chain direction is applied. The behavior of 
both setups is comparable as long as particles' displace- 
ments in x direction are small. However, this is no longer 
the case if the rearrangement to triangular structures has 
occurred. 

Finally, we used our simulation to track the time evo- 
lution of an aggregate comprising 200 primary particles 
in a resting fluid. The initial aggregate was obtained 
by diffusion- limited cluster particle aggregation [32| . We 
compared the results by using the classical DLVO forces 
only with the results obtained using our tangential force 
model. Figure [S] shows the initial aggregate and the re- 
structured aggregate after 25 seconds using DLVO forces 
only. It is remarkable that even in a resting fluid the 
aggregate collapses to more a compact structure. This 
behavior obviously is a result of the used force models. 
The van der Waals forces act over a relatively long range. 
As there is no resistance of single particle bonds against 




FIG. 6: Compaction of an aggregate using central forces only, 
(a) The start configuration, generated by diffusion limited 
aggregation, (b) The configuration after a simulated time 
period of 25s, the structure is significantly more compact than 
the starting configuration. 



bending moments, the bonded particles can freely re- 
orient. This behavior contradicts the observation that 
fractal structure are often stable over long times even in 
moderate shear flows [33, 34]. This indicates that assum- 
ing purely central forces is an over-simplification which 
hinders prediction of realistic structure evolution. Under 
the same conditions, but using the introduced tangential 
force model, the compaction does not occur. In order to 
quantify the compaction of the aggregate we count the 
number of neighboring particles for each particle. The 
average number of neighbors is used as a measure for the 
compactness. Particles are considered to be neighbors 
if their surface-surface distance is below 10 nm. Note 
that we deliberately do not use the fractal dimension as 
a measure for compactness because the aggregate looses 
the property of self-similarity during the compaction pro- 
cess. Figure [7] shows the time evolution of the average 
number of neighbors for the case of using DLVO forces 
only and for the case of using DLVO augmented with 
our tangential model. In case of using DLVO only, the 
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FIG. 7: Time evolution of the average number of neighbors. 
The solid line shows simulation results where DLVO forces 
only were considered. The dashed line displays data obtained 
from simulations where the tangential force model was ap- 
plied. 




FIG. 8: Sketch for the calculation of the bending energy. The 
initial configuration is shown by solid lines. After the de- 
formation process the three particles have new positions and 
orientations, shown by the dashed lines. The particles I and 
j are reoriented relative to particle i by the angles 8(j>i and 
8(f)i, and particle i was rotated around its center by the angle 
Sifi. The springs £ij and £il are deformed due to the particle 
motion. 



number of neighbors grows continuously, while the com- 
paction of the aggregate is negligible if the tangential 
model is used. 



APPENDIX A: DETERMINATION OF THE 
EQUILIBRIUM BENDING ENERGY IN LINEAR 
STRUCTURES 



VI. SUMMARY AND CONCLUSION 

In this work we have introduced a new model for tan- 
gential interaction forces applicable in computer simu- 
lations. Our tangential force model is capable of sup- 
porting bending moments. This is an important qual- 
ity as it has been recently shown experimentally that 
colloidal bonds actually support bending moments [2lj . 
Our model is based on two tangential springs acting be- 
tween bonded particles. The time evolution of these 
springs is determined by the relative tangential veloc- 
ity between the particles. The spring elongation is con- 
strained to a maximum elongation to account for sliding 
effects. The model contains two parameters: The stiff- 
ness of the springs and the maximum elongation. Wc 
showed that both parameters can be determined directly 
from the experiments shown in [2l| and simulations us- 
ing our model are able to reproduce the experimental 
observations. We compared our model to former models 
which were used in computer simulations and are capa- 
ble of supporting bending moments. Furthermore, we 
showed that even in the case of simulating a fractal ag- 
gregate in a resting fluid models without tangential forces 
lead to an unrealistic collapse of the aggregate. Using our 
tangential force model can remedy that flaw. 



In order to determine the bending energy in a linear 
structure of particles we consider three particles I, i and j 
of a structure which undergoes deformations as sketched 
in figure [5J The torque balance for particle i follows from 
equation (fl3|) 

= 2ak t x tin + £^ x ray) . (Al) 

Assuming small deformations, equation (|A1|) can be 
rewritten as a scalar equation: 

& - & = . (A2) 

Under the assumption that all deflections are small 
and no plastic deformation occurrs during deformation 
(i. e. no spring length in the system exceeded ^ ma x) one 
finds from equation (|13p : 

ta = 2a(6<tn-6<pi) , (A3) 
&j = 2a (5(f) j + Sifi) . (A4) 

Introducing this in equation (|A2[) and solving the result- 
ing equation for 6ipi leads to 

Sifi = hsfa -Sfa) . (A5) 

The energy stored in the springs £y and £y is given by 

Ei = \k t {^+^) (A6) 
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FIG. 9: Sketch for the calculation of the equilibrium shape 



Introducing (|A3|) and (|A4|) in (|A6|) leads to the equation 

Ei = 2a 2 h (6$ + 5<t> 2 + 28(p i (S<l> j - Sfa) + 28<p 2 ) . 

(A7) 



By introducing equation (|A5|) in (|A7|) one obtain: 

Ei = a 2 k t {5<j>i + <%) 2 . (A8) 

This in turn can be expressed as the change of the angle 
between the bond of particles I and i and particles i and 
j: 5$ = + 5(f>2- Therefore we get 

E, = a 2 k t 8<S> 2 . (A9) 

By adding up Ei for all particles in the chain one gets the 
whole elastic energy of the structure, except the contri- 
bution from rods rigidly connected to terminal particles 
to a particle within the structure. Therefore Ei can be 
understood as the energy per particle in an elastic de- 
formed linear structure of colloidal particles. 



APPENDIX B: DETERMINATION OF THE 
EQUILIBRIUM SHAPE 

In order to calculate the shape of a bended chain of 
particles as described in IIVB1 we minimize the elastic 
energy stored in such an elastic chain. From symmetry 
reasons it is sufficient to consider only half of the chain. 
As shown in figure [HI the change in the angle of the bonds 
connected to n'th particle is denoted as <f> n . The bending 
energy per particle is therefore given by (|A9[) : 

E n = a 2 k t 4> 2 n (Bl) 

The elastic energy stored in the whole chain is therefore 

JV-l 

E = 2 J2 a " k t€ » (B2) 



where N is the number of particles in the half of the 
chain. 

In equilibrium this energy should be minimal. The 
position of the n'th particle is 



Un 



-2a ^ sin ^ ^ 



i=l \j=i 
Assuming small total deflection this reduces to 



n i 



(B3) 



(B4) 



Assuming that the total deflection of the rod is As we 
have a constraint for the energy minimization: 



N-l i 
t=l j=l 



(B5) 



Using the technique of Lagrangian multiplier we have to 
minimize the quantity. 

N-l I N-i i \ 

£ = ]T 2a 2 fc t 02 _ A ha^J^+Aa (B6) 

i=l \ 1=1 j=l ) 

where A is the Lagrangian multiplier. To minimize E one 
has to solve 



dE 



N-l i 

= 4a 2 k t ^ m - 2Xa V V 5 3 . m , (B7) 

UVm t =l 3 = 1 

where Sj lTn is the Kronecker symbol which is equal to 
one if j = m and zero elsewhere. The solution of this 
equation is 



X(N -m) 
2ak t 



(B8) 



Introducing this in the constraint (|B5|) one finds after 
some algebra 



A fljv 3 - -N 2 + -N j : A.s 
h V3 2 6 ' 



(B9) 



Under neglecting all powers of N smaller then three one 
finds 



A = 



3fc t As 
N 3 



and by introducing that in equation (IB8P 

3 As ,„ 

(N - m) . 



2aN 3 



(BIO) 



(Bll) 



Introducing this in equation (|B4[) . using N = L/Aa, m = 
x m /a and neglecting all terms that vanish for a — > one 
finds 



Urn — ^3 (6x m 4\x m \ ) 



(B12) 



Introducing (|B8p in (|B2|) results in an expression for the 
elastic energy as function of N or L, respectively, and 
As. The bending force is related to that energy by 



dE 
dAs 



(B13) 



Now one can eliminate As. After introducing As in 
(|B12[) the elongation is found: 



-Pkend I L 2 1 i 



(B14) 



= 1 3=1 
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